C     MUNIMAT    Main subroutines.
C     Copyright (C) 1997  Filip Hroch, Masaryk University, Brno, CZ
!
!  $Id$
!                    
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

      FUNCTION UHEL( A,B )
C Arc between vectors. 
      
      DOUBLE PRECISION UHEL,B(2),A(2),C2,A2,B2,X
      
      C2 = ( A(1)*B(2)-A(2)*B(1) )**2
      A2 = A(1)**2+A(2)**2
      B2 = B(1)**2+B(2)**2                   
      IF( A2*B2.EQ.0.0 )THEN
        UHEL = -1.0
      ELSE
        X = SQRT(C2/(A2*B2))
        IF( ABS(X).GT.1.0 ) X = 1.0
        UHEL = ASIN(X)
      ENDIF
      
      END
      
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      
      LOGICAL FUNCTION MATCH (N,ID1,ID2,X1,Y1,X2,Y2,C)  
C Do is correlation between a stars in series ID1 and ID2 ?     

      INCLUDE 'munimat.inc'                     
      
      INTEGER IDMAX2
      DOUBLE PRECISION XTOL
      PARAMETER( IDMAX2 = IDMAX*IDMAX/2+1, XTOL = 0.05 )
      
      INTEGER N,I1,I2,I3,ID1(*),ID2(*),I,J,K,L,IPOM(IDMAX2)
      DOUBLE PRECISION 
     . X1(*),Y1(*),X2(*),Y2(*),K2(IDMAX2),R21,R22,XMED,S,XDEV
      DOUBLE PRECISION A1(2),A2(2),B1(2),B2(2),UHEL
      DOUBLE PRECISION PRUM,SIG,DEV(IDMAX2),XARC,MAXDEV,C
      LOGICAL REP

      MATCH = .FALSE.          
C Tests on the same length sides of N-polygon.
C Calculation of ratio of similarity (scale).
      REP = .FALSE.
      L = 0
      DO I=1,N
        DO J=I+1,N  
          L = L+1
          I1 = ID1(I)
          I2 = ID1(J)
          R21 = (X1(I2) - X1(I1))**2 + (Y1(I2) - Y1(I1))**2
          I1 = ID2(I)
          I2 = ID2(J)
          R22 = (X2(I2)-X2(I1))**2 + (Y2(I2)-Y2(I1))**2
          K2(L) = SQRT(R22/R21)
          IF( L.EQ.1.AND.K2(L).LT.0.9 ) REP = .TRUE. 
          IF( REP ) K2(L) = 1.0/K2(L)
        ENDDO  
      ENDDO
      CALL USPOR(L,K2,IPOM)
C Calculation of medians.
      IF( MOD(L,2).EQ.0 )THEN
        XMED = 0.5*(K2(L/2)+K2(L/2+1))
      ELSE
        XMED = K2(L/2+1)
      ENDIF                    
c      print *,'meritka=',(k2(i),i=1,l)  
c      print *,'id1>',(id1(i),i=1,n)
c      print *,'id2>',(id2(i),i=1,n)

C Calculation of absolute deviations.
      DO I=1,L
        DEV(I) = ABS(K2(I)-XMED)
C  a deviation is too great (because it isn't from the true series)
        IF( DEV(I)/XMED.GT.0.2 ) RETURN
      ENDDO
      CALL USPOR(L,DEV,IPOM)
C Calculation of the median deviation => calculation of the mean deviation
      IF( MOD(L,2).EQ.0 )THEN
        XDEV = 0.5*(DEV(L/2)+DEV(L/2+1))
      ELSE
        XDEV = DEV(L/2+1)
      ENDIF                      
C Mean deviation is too great
c      print *,xdev/xmed
      IF( XDEV/XMED.GT.XTOL ) RETURN
c        print *,xmed,k2(1),k2(l)
C Calculation of the aritmetical mean
      S = 0.0
      DO I=1,L
        S = S+K2(I)        
      ENDDO
      PRUM = S/L  
C Find the greatest deviation into median
      J = 0
      MAXDEV = -1.0
      DO I=1,L
        DEV(I) = (K2(I)-PRUM)**2
        IF( DEV(I).GT.MAXDEV )THEN
          MAXDEV = DEV(I)
          J = I   
        ENDIF
      ENDDO        
C Calculation of the new mean.
      S = 0.0
      DO I=1,L
        IF( I.NE.J ) S = S+K2(I)        
      ENDDO
      PRUM = S/(L-1)  
C Estimation of the standart deviation
      S = 0.0
      DO I=1,L
        DEV(I) = K2(I) - PRUM
        IF( I.NE.J ) S = S+DEV(I)**2
      ENDDO
      SIG = SQRT(S/(L-2))+MACHEPS
C Standard deviation is too great.
      IF( SIG/PRUM.GT.XTOL ) RETURN
C Test for equality of median and mean
      IF( ABS(PRUM-XMED).GE.SIG ) RETURN          

c      PRINT *,XMED,XDEV,PRUM,SIG          

C Test on the same arcs of N-polygon
      L = 0
      DO I=1,N
        DO J=I+1,N
          DO K=J+1,N
            L = L+1
            I1 = ID1(I)
            I2 = ID1(J)
            I3 = ID1(K)
            A1(1) = X1(I2)-X1(I1)
            A1(2) = Y1(I2)-Y1(I1)
            B1(1) = X1(I3)-X1(I1)
            B1(2) = Y1(I3)-Y1(I1)
            I1 = ID2(I)
            I2 = ID2(J)
            I3 = ID2(K)
            A2(1) = X2(I2)-X2(I1)
            A2(2) = Y2(I2)-Y2(I1)
            B2(1) = X2(I3)-X2(I1)
            B2(2) = Y2(I3)-Y2(I1)
            K2(L) = UHEL(A1,B1)/UHEL(A2,B2)
C          devation onto 1.0 is too great            
c          print *,'.',k2(l)
            IF( K2(L).GT.2.0 .OR. K2(L).LT.0.5 ) RETURN
          ENDDO
        ENDDO
      ENDDO      
      CALL USPOR(L,K2,IPOM)
C Calculation of the median of arcs
      IF( MOD(L,2).EQ.0 )THEN
        XARC = 0.5*(K2(L/2)+K2(L/2+1))
      ELSE
        XARC = K2(L/2+1)
      ENDIF                      
c      print *,xarc
C Calculation of the deviation into median
      DO I=1,L
        DEV(I) = ABS(K2(I) - XARC)
      ENDDO
      CALL USPOR(L,DEV,IPOM)
C Calculation of the median from deviation 
      IF( MOD(L,2).EQ.0 )THEN
        XDEV = 0.5*(DEV(L/2)+DEV(L/2+1))
      ELSE
        XDEV = DEV(L/2+1)
      ENDIF                      
      XDEV = XDEV + MACHEPS
C Mean deviation is too great
c     print *,'..',xdev,xarc
      IF( XDEV/XARC.GT.0.2 ) RETURN
c      print *,xdev,xarc
C Median must be equal 1.0
      IF( ABS(XARC-1.0).GT.3.0*XDEV ) RETURN
c      print *,xarc,xdev,k2(1),k2(l)

c     print *,xarc,xdev

C All test are OK!        
      MATCH = .TRUE.
      IF( REP )THEN
        C = 1.0/XMED
      ELSE   
        C = XMED
      ENDIF  
c      print *,sqrt(c),sig1
      END

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

      SUBROUTINE INIT1(I,N1,X,Y,N,ID,I1,I2,Z1)

      INCLUDE 'munimat.inc'
      INTEGER N1,N,ID(*),I,J,I1,I2,IDR(MMAX)
      DOUBLE PRECISION X(*),Y(*),Z1 
      DOUBLE PRECISION XC,YC,R2(MMAX)
      COMMON /IXX1/IDR,R2
      SAVE /IXX1/
                                                
C Initialization     
      XC = 0.5*(X(I1) + X(I2))
      YC = 0.5*(Y(I1) + Y(I2))
      DO J=1,N1                    
        R2(J) = (X(J) - XC)**2+(Y(J) - YC)**2
        IDR(J) = J
        IF( J.EQ.I1 ) R2(J) = -2.0
        IF( J.EQ.I2 ) R2(J) = -1.0
      ENDDO
      CALL USPOR2(N1,R2,IDR)                                         
      ID(1) = IDR(1)
      ID(2) = IDR(2)
      Z1 = (X(I1) - X(I2))**2 + (Y(I1) - Y(I2))**2
      N = 2              
      I = 2
      END
              
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

      SUBROUTINE INIT2(I,N1,X,Y,N,ID,I1,I2,Z1)

      INCLUDE 'munimat.inc'
      INTEGER N1,N,ID(*),I,J,I1,I2,IDR(MMAX)
      DOUBLE PRECISION X(*),Y(*),Z1
      DOUBLE PRECISION XC,YC,ZOOM,Z2,R2(MMAX)
      COMMON /IXX2/IDR,R2
      SAVE /IXX2/                                          

C Initialization
      XC = 0.5*(X(I1) + X(I2))
      YC = 0.5*(Y(I1) + Y(I2))
      Z2 = (X(I1) - X(I2))**2 + (Y(I1) - Y(I2))**2     
      ZOOM = Z1/Z2
c      RT = 2.0*4.0*TOL**2*ZOOM*(1.0/Z2+Z1/Z2**2)
c      RT = 4.0*16.0*TOL**2
      DO J=1,N1                    
        R2(J) = ZOOM*((X(J) - XC)**2 + (Y(J) - YC)**2)
        IDR(J) = J            
        IF( J.EQ.I1 ) R2(J) = -2.0
        IF( J.EQ.I2 ) R2(J) = -1.0
      ENDDO
      CALL USPOR2(N1,R2,IDR)                                         
      ID(1) = IDR(1)
      ID(2) = IDR(2)
      N = 2              
      I = 3
      END

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

      SUBROUTINE UV(N,ID,X,Y,U,V)
      
      INTEGER N,ID(*),L,J,K
      DOUBLE PRECISION X(*),Y(*),U,V
      DOUBLE PRECISION D(3)
      
      L = ID(N-2)
      J = ID(N-1)
      K = ID(N)                        
      D(1) = (X(L)-X(J))**2+(Y(L)-Y(J))**2
      D(2) = (X(L)-X(K))**2+(Y(L)-Y(K))**2
      D(3) = (X(K)-X(J))**2+(Y(K)-Y(J))**2
      CALL SERAD(D)
      U = SQRT(D(2)/D(1))
      V = SQRT(D(3)/D(1))
      END

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

      LOGICAL FUNCTION FIND1(I,N1,N,ID1,R)
C Find next candidate
      INCLUDE 'munimat.inc'
      INTEGER I,N,ID1(*),IDR(MMAX),N1
      DOUBLE PRECISION R2(MMAX),R
      COMMON /IXX1/IDR,R2
      SAVE /IXX1/                                          
                                         
      IF( I.GE.N1.OR.N.GE.IDMAX ) THEN                                     
        FIND1 = .FALSE.                                                      
      ELSE  
        N = N+1   
        I = I+1
        ID1(N) = IDR(I)                
        R = R2(I)
        FIND1 = .TRUE.
      ENDIF  
      END      
      
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC                                      

      LOGICAL FUNCTION FIND2(I,N2,X2,Y2,IN2,ID2,R,N1,X1,Y1,IN1,ID1)

      INCLUDE 'munimat.inc'
      INTEGER IN1,IN2,I,I0,J,L,K,ID2(*),ID1(*),N1,N2,IDR(MMAX)
      DOUBLE PRECISION X1(*),Y1(*),X2(*),Y2(*),RL,RG,D(3)
      DOUBLE PRECISION U1,V1,U2,V2
      DOUBLE PRECISION R,R2(MMAX),RMAX,DIF
      LOGICAL PROJIT 
      COMMON /IXX2/IDR,R2
      SAVE /IXX2/                                          
      
      RL = 0.7*R              
      RG = 1.3*R    
      PROJIT = .TRUE.

C Find first star for the condition  RL < R2(I) < RG
5     IF( R2(I-1).GE.RL.AND.IDR(I-1).NE.ID2(IN2)
     .   .AND.I.LT.N2.AND.I.GT.3 )THEN
        I = I-1
        GOTO 5
      ENDIF  

C Find star in the annulus      
      I0 = 0                           
      RMAX = MACHMAX                           
10    IF( I.LE.N2.AND.R2(I).LE.RG )THEN
        DIF = ABS(R2(I)-R)
        IF( DIF.LT.RMAX )THEN
          I0 = I
          RMAX = DIF
        ENDIF     
        I = I+1
        GOTO 10 
      ENDIF     
      
      IF( I0.EQ.0 )THEN
        FIND2 = .FALSE.                                             
      ELSE  
        FIND2 = .TRUE.                                             
        I = I0 
C Add star to series
        IN2 = IN2+1
        ID2(IN2) = IDR(I)             
C Calculation of the triangle in the first file        
        CALL UV(IN1,ID1,X1,Y1,U1,V1)   
C Calculation of the triangle in the second file        
        L = ID2(IN2-2)
        J = ID2(IN2-1)
        K = ID2(IN2)
20      D(1) = (X2(L)-X2(J))**2+(Y2(L)-Y2(J))**2
        D(2) = (X2(L)-X2(K))**2+(Y2(L)-Y2(K))**2
        D(3) = (X2(K)-X2(J))**2+(Y2(K)-Y2(J))**2
        CALL SERAD(D)
        U2 = SQRT(D(2)/D(1))
        V2 = SQRT(D(3)/D(1))
        IF( (U1-U2)**2+(V1-V2)**2.LT.5E-3 )THEN
C Transformation is found.
          I = I+1    
c         if( id1(1).eq.2.and.id1(2).eq.3)then
c           print *,'u:',u1,u2,abs(u1-u2)
c           print *,'v:',v1,v2,abs(v1-v2)
c         endif
          IF( .NOT.PROJIT )THEN
            L = ID2(IN2-2)
            ID2(IN2-2) = ID2(IN2-3)
            ID2(IN2-3) = L
          ENDIF  
          RETURN
        ENDIF
C Case, that firts two stars in triangle are turn between files   
        IF( IN2.EQ.4.AND.PROJIT )THEN
          L = ID2(IN2-3)
          PROJIT = .FALSE.
          GOTO 20
        ENDIF  
C star is in annulus, but U,V are a bad.
        IN2 = IN2-1
        FIND2 = .FALSE.
      ENDIF
      END

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

      SUBROUTINE SOLVE(WATCH,NSTAR,CLIP,N1,X1,Y1,N2,X2,Y2,IP,S,R,M,C,N)
     
      INCLUDE 'munimat.inc'
      DOUBLE PRECISION MACH9                    
      INTEGER NSTAR
      PARAMETER( MACH9 = 0.9*MACHMAX )
                      
      INTEGER I,J        
      INTEGER N,N1,N2
      DOUBLE PRECISION 
     .  X1(*),Y1(*),X2(*),Y2(*),S,R(2),M(2,2),C,Z1,TRAFO,SMIN
      DOUBLE PRECISION SSTACK(MAXSTACK),CLIP
      DOUBLE PRECISION RR
      LOGICAL IP
      INTEGER LL,NI1,NI2,ID1(IDMAX),ID2(IDMAX),II1,II2
      INTEGER I1,I2,J1,J2
      INTEGER NSTACK,STACK1(MAXSTACK,0:IDMAX),STACK2(MAXSTACK,0:IDMAX)
      LOGICAL FIND1,FIND2,MATCH,WATCH
                          
      NSTACK = 0                            
      SMIN = MACHMAX
      DO I1=1,N1       
        IF( WATCH ) PRINT *,'INDEX = ',I1
        DO 95,I2=I1+1,N1                
          CALL INIT1(II1,N1,X1,Y1,NI1,ID1,I1,I2,Z1)
          DO J1=1,N2
            DO 90,J2=J1+1,N2 
              II1 = 2                     
              NI1 = 2
              CALL INIT2(II2,N2,X2,Y2,NI2,ID2,J1,J2,Z1)
80            IF( FIND1(II1,N1,NI1,ID1,RR) )THEN
                IF(FIND2(II2,N2,X2,Y2,NI2,ID2,RR,N1,X1,Y1,NI1,ID1))THEN
                  IF( NI1.GE.NSTAR.AND.NI2.GE.NSTAR )THEN 
                    IF( .NOT.MATCH(NSTAR,ID1,ID2,X1,Y1,X2,Y2,C) )
     *                GOTO 90      
                    S = TRAFO(NSTAR,CLIP,
     .                        ID1,ID2,N1,N2,X1,Y1,X2,Y2,R,M,C,N)
                    IF( S.LT.1E3*SMIN )THEN
C Transformation is find.
C They parameters is best as previously, recording the transformation
                      IF( WATCH )THEN 
                        PRINT 2,' FIN1>',(ID1(LL),LL=1,NI1)
                        PRINT 2,' FIN2>',(ID2(LL),LL=1,NI2)
                        PRINT 1,S,N,R(1),R(2),M(1,1),M(1,2),C
                      ENDIF
                      IP = .TRUE.
                      SMIN = MIN(SMIN,S)
                      CALL ADD(NSTAR,ID1,ID2,NSTACK,STACK1,
     .                         STACK2,SSTACK,S)
                    ENDIF
                    GOTO 90
                  ENDIF
                ELSE
                  NI1 = NI1-1
                  IF( II2.GT.N2 )THEN
                    NI1 = 2
                    GOTO 90
                  ENDIF  
                ENDIF
                GOTO 80
              ENDIF         
90          CONTINUE  
          ENDDO  
95      CONTINUE
      ENDDO   

      IF( .NOT.IP ) RETURN

      IF( WATCH )THEN
        PRINT *,' Nstack=',nstack
        DO J=1,NSTACK
          IF( SSTACK(J).LT.1E3*SMIN )THEN
            PRINT *,'1>',(STACK1(J,I),I=0,NSTAR)
            PRINT *,'2>',(STACK2(J,I),I=0,NSTAR)
          ENDIF
        ENDDO  
      ENDIF
C Find the most numerous series and calculation of result transformation
      J = 0  
      DO I=1,NSTACK
        IF( SSTACK(I).LT.1E3*SMIN )THEN
          IF( STACK1(I,0).GT.J )THEN
            J = STACK1(I,0)
            LL = I
          ENDIF
        ENDIF
      ENDDO  
      IF( J.LE.2 ) THEN
        PRINT *,'I''m a somewhat confused...sorry,',
     *          ' but this I can''t identify.'
        IP = .FALSE.
        RETURN
      ENDIF
      DO I=1,NSTAR
        ID1(I) = STACK1(LL,I)
        ID2(I) = STACK2(LL,I)
      ENDDO
      IF( WATCH )THEN  
        PRINT *,' Cross references in the most numerous case : ',
     *        STACK1(LL,0)
        PRINT 2,' FIN1>',(ID1(LL),LL=1,NSTAR)
        PRINT 2,' FIN2>',(ID2(LL),LL=1,NSTAR)
      ENDIF
      IF( .NOT.MATCH(NSTAR,ID1,ID2,X1,Y1,X2,Y2,C) ) STOP 'Error 010'
      S = TRAFO(NSTAR,CLIP,ID1,ID2,N1,N2,X1,Y1,X2,Y2,R,M,C,N) 
      IF( S.GE.MACH9 ) STOP 'Error 020'
      
1     FORMAT(1P,E10.3,0P,I4.1,2F8.2,2F9.5,1P,E12.5)
2     FORMAT(A,90I5)
      END

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

      SUBROUTINE TRI(GETARG1,GETARG2,PARY,trapar)

      INCLUDE 'munimat.inc'

      DOUBLE PRECISION RAD
      PARAMETER( RAD = 180.0/3.14159 26535 89793 23846 )

      DOUBLE PRECISION 
     .  X1(MMAX),Y1(MMAX),X2(MMAX),Y2(MMAX),R(2),M(2,2)
      DOUBLE PRECISION S,C,T1(2),T2(2),X
      DOUBLE PRECISION CLIP
      DOUBLE PRECISION TOL,RMIN,DX,DY,DR,XX,YY
      CHARACTER*80 SOUBOR1,SOUBOR2,GETARG1,GETARG2,SWITCH*250
      character*80 fsoubor1, fsoubor2
      CHARACTER*10 PARY,trapar
      CHARACTER*256 RADEK1,RADEK2
      INTEGER I,J,J0,N,N1,N2,IN1(MMAX),IN2(MMAX),II,L,LN
      INTEGER MAXSTAR,NSTAR,NL
      LOGICAL IP,WATCH,BLANKA,WRIBL,KOMOUT,CAT
      integer lineid1(MMAX), lineid2(MMAX)
      character*256 lines2(3*MMAX+3)

      CAT = PARY(1:1) .EQ. 'C'
      KOMOUT = PARY(2:2) .EQ. 'K'
C Read parameters from the file "match.opt"
      CALL GETPAR(MAXSTAR,NSTAR,CLIP,WATCH,BLANKA)           
      
      CALL CTOF(GETARG1,SOUBOR1)
      CALL CTOF(GETARG2,SOUBOR2)

      I = LN(SOUBOR1)
      J = LN(SOUBOR2)
      PRINT *,'Match files : ',SOUBOR1(:I),' -> ',SOUBOR2(:J)      
        
C Read data from SOUBOR1 and SOUBOR2
      CALL DATA(WATCH,SOUBOR1,N1,X1,Y1,IN1,MAXSTAR)               
      CALL DATA(WATCH,SOUBOR2,N2,X2,Y2,IN2,MAXSTAR)               
      IF( N1.LT.IDMAX.OR.N2.LT.IDMAX )THEN
        PRINT *,' Very few points.'
        RETURN
      ENDIF   
      IF( CAT )THEN
        DO I = 1,N1
          X1(I) = X1(I)/RAD
          Y1(I) = Y1(I)/RAD
        ENDDO
      ENDIF 
C Calcualtion of the mean of the coordinates 
      CALL CENTER(N1,X1,Y1,T1)
      CALL CENTER(N2,X2,Y2,T2)
      IF( CAT )THEN
        DO I = 1,N1
          X1(I) = (X1(i) - T1(1))*COS(t1(2))
          Y1(I) = Y1(I) - T1(2)
        ENDDO
      ELSE
        DO I = 1,N2
         X1(I) = X1(I) - T1(1)
         Y1(I) = Y1(I) - T1(2)
        ENDDO
      ENDIF
C pre-calculation coordinates of the second file 
      DO I = 1,N2
       X2(I) = X2(I) - T2(1)
       Y2(I) = Y2(I) - T2(2)
      ENDDO

      IF( WATCH )
     *PRINT *,'     SUM   PR      DX      DY',
     *    '   M(1,1)   M(1,2)   Scale'
      IP = .FALSE.
      
      CALL SOLVE(WATCH,NSTAR,CLIP,N1,X1,Y1,N2,X2,Y2,IP,S,R,M,C,N )

      IF( .NOT.IP )THEN
        PRINT *,' ******   Coincidences not found !  ******'
        RETURN
      ENDIF
      IF( WATCH )THEN
        PRINT '(A,I10.1)',' Number of matches stars : ',N
        PRINT '(A,1P,G20.5)',' Sum of square in the betters case : ',
     -                     S*(N-4)
        PRINT '(A,2F15.7)',' Mean 1: ',(T1(I),I=1,2) 
        PRINT '(A,2F15.7)',' Mean 2: ',(T2(I),I=1,2)
      ENDIF
      PRINT *
      PRINT '(A,2F15.10)','   Matrix : ',(M(1,I),I=1,2)
      PRINT '(A,2F15.10)','            ',(M(2,I),I=1,2)
      PRINT '(A,2F30.10)','   Shift  : ',(R(I),I=1,2)
      IF( C.GE.0.99 )THEN 
        PRINT '(A,1P,G20.10)','   Scale  : ',C
      ELSE
        PRINT '(2(A,1P,G20.10))','   Scale : ',C,' = 1/',1.0/C
      ENDIF

c vypsani parametru do souboru
      if( trapar .ne. ' ' )then
         open(31,file=switch(soubor2,'.TRA'),iostat=i)
         if( i .ne. 0 ) stop '31'
         fsoubor1 = switch(soubor1,'')
         fsoubor2 = switch(soubor2,'')
         i = ln(fsoubor1)
         j = ln(fsoubor2)
         write(31,'(5A,12G20.10)') '''',fsoubor1(:i),''' ''',
     -        fsoubor2(:j),'''',c, r, m, t1, t2
         close(31)
      endif

                         
C CTENI DAT ZE SOUBORU                     
      CALL DATA(WATCH,SOUBOR1,N1,X1,Y1,IN1,MMAX)               
      CALL DATA(WATCH,SOUBOR2,N2,X2,Y2,IN2,MMAX)               

! index the stars in soubor1 (the 'lineid1' contains line number of x1,y1)
      open(1,file=soubor1,status='old',blank='zero')
      i = 0
      j = 1
      do while(.true.)
         read(1,'(A)',end=970,err=970) RADEK1
         read(RADEK1,'(T7,2F9.3)',iostat=ii) x,y
         i = i + 1
         if( ii .eq. 0 ) then
            if( (x - x1(j))**2 + (y - y1(j))**2 .lt. 1e-4 ) then
               lineid1(j) = i
               j = j + 1
            endif
         endif
      enddo
 970  continue
      close(1)
      do i = 1, MMAX
         lineid2(i) = 0
      enddo

      TOL = S*CLIP**2
      IF( WATCH ) PRINT '(A,1P,G15.3)',' Tolerance :  +-',SQRT(TOL)
      PRINT '(1H ,3A)','      Id1        X1        Y1',
     *                 '       Id2        X2        Y2',
     *                 '       dx       dy' 
      DO I=1,N1
        IF( CAT )THEN
          XX1 = (X1(I)/RAD - T1(1))*COS(T1(2))
          YY1 =  Y1(I)/RAD - T1(2)
        ELSE
          XX1 = X1(I) - T1(1)
          YY1 = Y1(I) - T1(2)
        ENDIF
        XX = C*(R(1) + XX1*M(1,1) + Yy1*M(1,2))
        YY = C*(R(2) + XX1*M(2,1) + Yy1*M(2,2))
        RMIN = MACHMAX
        J0 = 1
        DO J=1,N2
          DX = XX - (X2(J) - T2(1))
          DY = YY - (Y2(J) - T2(2))
          DR = DX**2+DY**2
          IF( DR.LT.RMIN ) THEN
            RMIN = DR
            J0 = J
          ENDIF
        ENDDO
        J = J0
        DX = XX - (X2(J) - T2(1))
        DY = YY - (Y2(J) - T2(2))
        IF( RMIN.LT.TOL )THEN
          PRINT '(2(I10.1,2F10.3),2F9.5)',
     *      IN1(I),X1(I),Y1(I),IN2(J),X2(J),Y2(J),DX,DY
          lineid2(j) = i
       endif
      ENDDO

! make output .MAT file
      OPEN(1,FILE=SOUBOR2,STATUS='OLD',BLANK='ZERO',IOSTAT=I)
      IF( I.NE.0 ) STOP '050'                   
      SOUBOR2 = SWITCH(SOUBOR2,'.MAT')
      OPEN(2,FILE=SOUBOR2,IOSTAT=I)
      IF( I.NE.0 ) STOP 'Fatal error during open the MATch file.'
      IF( .NOT.KOMOUT )THEN
        CALL RDHEAD(1, NL, I, I, X, X, X, X, X, X, X)
        IF( NL.NE.0 )THEN
          REWIND(1)
          DO I=1,3
            READ(1,'(A)',IOSTAT=II) RADEK1
            IF( II.NE.0 ) STOP 'Fatal error for reading header of file.'
            WRITE(2,'(A)',IOSTAT=II) RADEK1(:LN(RADEK1))
            IF( II.NE.0 ) STOP 'Fatal error for reading header of file.'
          ENDDO
        ENDIF
      ENDIF
      WRIBL = .TRUE.

! header
      REWIND(1)
      IF( NL.EQ.0 )THEN
         J0 = 0
      ELSE
         J0 = 3
      ENDIF
      DO L=1,J0
         READ(1,*,IOSTAT=II) 
         IF( II.NE.0 ) STOP 'Fatal error during read file1.'
      ENDDO
      READ(1,'(A)',IOSTAT=II) RADEK1
      IF( II.NE.0 ) STOP 'Fatal error during read file2.'
      IF( RADEK1.NE.' ' )THEN
         BACKSPACE(1,IOSTAT=II)
      ELSE
         IF( WRIBL ) WRITE(2,*,IOSTAT=II)
      ENDIF
      WRIBL = .FALSE.
      IF( II.NE.0 ) STOP 'Fatal error during read file3.'

      do i = 1, 3*MMAX
         lines2(i) = ' ' 
      enddo

! file to match read
      maxir = 0
      do j = 1, n2
! reindex stars
         if( 1 .le. lineid2(j) .and. lineid2(j) .le. n1 ) then
         ir = lineid1(lineid2(j))
         IF( NL.EQ.2 )THEN
            READ(1,'(A/A/)',END=95,IOSTAT=II) lines2(ir),lines2(ir+1)
         ELSE
            READ(1,'(A)',END=95,IOSTAT=II) lines2(ir)
         ENDIF

! output modes
         IF( KOMOUT )THEN
            WRITE(lines2(ir),'(I6,2F9.3)',IOSTAT=II) IN1(I),X2(J),Y2(J)
         ENDIF      
         IF( II.NE.0 )  STOP 'Fatal error during read file4.'   
         if( ir .gt. maxir ) maxir = ir
         endif
      enddo
 95   continue

      do j = lineid1(1),maxir+2
         write(2,'(A)') lines2(j)(:ln(lines2(j)))
      enddo

      CLOSE(1)
      CLOSE(2)
      
      END

